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ABSTRACT 

, A new method, based on a Reverse Monte Carlo technique and aimed at the inverse 

' problem in the analysis of interstellar (intergalactic) absorption lines is presented. 

We consider the process of line formation in media with a stochastic velocity field 
t-H ■ accounting for the effects of a finite correlation length (mesoturbulence) . This approach 

generalizes the standard microturbulent approximation, which is commonly used to 
s ! , model the formation of absorption spectra in turbulent media. 

(N . The method allows to estimate from an observed spectrum both, the physical pa- 

| rameters of the absorbing gas, and an appropriate structure of the velocity distribution 

parallel to the line of sight. The procedure is applied to a template H+D Lya profile 
which reproduces the Buries & Tytler [B&T] spectrum of the quasar Q 1009+2956 
| with the DI Lya line seen at the redshift z a = 2.504. 

The results obtained favor a low D/H ratio in this particular absorption system, 
' although the inferred upper limit for the hydrogen isotopic ratio of about 4.6 x 10~ 5 is 

slightly higher than that of B&T [log (D/H) = -4.60 ± 0.08 ± 0.06, i.e. D/H^ 3.47 x 
Q ■ 10~ 5 ]. This revision of the B&T data in the framework of the mesoturbulent model 

\ leads to limits on D/H consistent with standard big bang nucleosynthesis [SBBN] 

c/3 . predictions and observational constraints on both extra-galactic 4 He mass fraction Y p 

(Izotov et al. 1997) and 7 Li abundance in the atmospheres of population II (halo) 
J> \ stars (Bonifacio & Molaro 1997). SBBN restricts the common interval of the baryon 

• ■ to photon ratio rj for which there is concordance between the various abundances to 

p>S ; 3.9 x lO" 10 ^i|^5.3x lO" 10 . It implies that 0.014 ^ Q, b h 2 <> 0.020, where h is the 

^ , Hubble constant in units 100 km s -1 Mpc -1 . 

Key words: line: formation - line: profiles - IGM: absorption lines - abundances: 
D/H - quasars: absorption lines - quasars: individual: Q 1009+2956. 



1 INTRODUCTION 

In two previous papers (Levshakov & Kegel 1997, and Lev- 
shakov, Kegel & Mazets 1997, Paper I and II hereinafter, 
respectively), we investigated the formation of absorption 
lines in chaotic media with finite correlation length of the 
velocity field [mesoturbulence), when the observed absorp- 
tion spectra correspond to only one line of sight. We have 
shown that accounting for a finite correlation length in the 
large scale stochastic velocity field changes the interpreta- 
tion of absorption spectra considerably. The line profiles, as 
well as the equivalent widths, may be strongly affected by 
the correlations. Therefore the standard analysis of absorp- 
tion spectra based on the microturbulent approximation and 



the Voigt profile fitting procedure may yield misleading re- 
sults concerning the physical parameters of the absorbing 
gas (for numerous examples, see Levshakov & Kegel 1994, 
1996; Levshakov & Takahara 1996a, b; Paper I and Paper II). 
At least one can say that the results obtained by the stan- 
dard analysis are model dependent and, thus, are not unique. 

In our cloud model we account for a stochastic velocity 
field with finite correlation length I but assume the absorb- 
ing material to be homogeneous in density and tempera- 
ture. The model is fully defined by specifying the column 
density N, the kinetic temperature Thin, the ratio of the 
r.m.s. turbulent velocity to the thermal velocity at/vth, the 
ratio of the cloud thickness to the correlation length L/l, 
as well as one realization of the velocity field distribution 
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along the line of sight. Because of the very large time scale 
for changes in the hydrodynamic flows the random structure 
of the velocity field along a given line of sight has to be con- 
sidered as 'frozen' over the exposure time. It follows that 
the observed absorption spectrum in the light of a point- 
like source (star, QSO) corresponds to only one realization 
of the velocity field. Therefore, the intervening absorption 
spectrum cannot be considered as a time or space average 
in the statistical sense. This means that the actual distri- 
bution of the hydrodynamic velocities at a given instant of 
time corresponds to an incomplete sample and, thus, may 
deviate from the average distribution function assumed to 
be Gaussian. Accordingly, significant deviations from the ex- 
pected average intensity (7a) can occur. - For this reason we 
used in Paper II a Monte Carlo technique to calculate spec- 
tra corresponding to the absorption along individual lines of 
sight. 

In Paper II we considered the direct problem, i.e. we 
specified the physical parameters and generated individual 
random realizations of the velocity distribution [more ex- 
actly - the distribution of the velocity component parallel to 
the line of sight v(s)] with which we then calculated individ- 
ual spectra. The aim of the present paper is the inverse prob- 
lem, i.e. the problem to deduce physical parameters from an 
observed spectrum. To reproduce a given (observed) spec- 
trum within the framework of our model, one has to find the 
proper physical parameters and an appropriate realization 
of the velocity field. In general, v(s) is a continuous random 
function of the coordinate s, but in the numerical procedure 
it is sampled at evenly spaced intervals As. The necessary 
number of intervals depends on the values of at/vth and L/l, 
being typically ~ 100 for hydrogen absorption lines. Thus, 
to estimate model parameters from the observed spectrum 
one has to solve an optimization problem in a parameter 
space of very large and variable (depending on a t /v t h, L/l) 
dimension. It is known that such kind of problems may be 
solved using stochastic optimization methods. The Reverse 
Monte Carlo [RMC] technique based on the computational 
scheme invented by Metropolis et al. (1953) proved to be ad- 
equate in our case. The method is successfully used in many 
applications where the state space of a physical system is 
huge (see e.g. Press et al. 1992, or Hoffmann 1995). 

We apply the RMC method to the problem of determin- 
ing the primordial deuterium abundance at high redshift. In 
particular, the H+D Lya profile observed by Buries & Tytler 
(1996, B&T hereinafter) in the spectrum of the quasar Q 
1009+2956 is considered. This spectrum was selected since 
(a) it shows a pronounced DI absorptions at the redshift 
z a = 2.504, (b) it was obtained with high signal-to-noise ra- 
tio and spectral resolution, and (c) the total hydrogen col- 
umn density estimated by B&T from the intensity level be- 
yond the Lyman limit in this system is consistent with the 
normalized intensities in the HI Lya wings. 

In the present paper we extend our study of the ac- 
curacy of the D/H ratio determinations from our Galaxy 
(Paper II) to the very distant Lya-systems (putative inter- 
vening galaxies). The background source is supposed to be 
point-like and even in the case of distant QSOs we consider 
absorption along one line of sight only. This means that 
any gravitational focusing which may, in principle, lead to 
an additional spatial averaging over different lines of sight 
(an example may be found in Frye et al. 1997) is not taken 



into account. (The average mesoturbulent H+D Lya spectra 
have been considered by Levshakov & Takahara 1996a.) 

In Section 2 a description of the RMC method and the 
procedure to test its validity are given, while in Section 3 
the method is applied to analyze an observed H+D Lya 
spectrum. The results obtained are summarized in Section 



2 THE REVERSE MONTE-CARLO METHOD 

The RMC method belongs to the class of stochastic opti- 
mization algorithms developed to solve optimization prob- 
lems with a very large number of free parameters. Contrary 
to the standard Monte-Carlo procedure in which random 
configurations of a given physical system are generated to 
estimate its average characteristics, the RMC takes an ex- 
perimentally determined set of data and searches for a pa- 
rameter configuration which reproduces the observational 
data. 

The inverse problem is always an optimization prob- 
lem in which an objective function is minimized. The objec- 
tive function defined over a space of very large dimension 
is known to have many local minima. Therefore optimiza- 
tion methods based on the choice of 'minimization direction' 
(so-called gradient methods) may fail because of the possi- 
bility of being trapped in one of such minima. The RMC 
approach allows to overcome this problem. It can get over 
the local minima because it does not use any concept of 'di- 
rection'. Applied to the analysis of absorption spectra, the 
RMC method may be formulated as follows. 



2.1 Computational procedure 

The total parameter space consisting of the physical parame- 
ters (like N, Thin, etc.) and the velocity components Vj par- 
allel to the line of sight at the individual positions Sj, is 
divided into two subspaces. 

Let O — {9} represent the subspace of physical parame- 
ters {-/V, Tkin, o~t/vth, L/l} . For the special case of H+D Lya 
that we want to investigate, 9 has to have one additional 
component D/H - the ratio of the DI to HI column densi- 
ties. 

Let further {v} = {vi,V2, ■ ■ ■ ,Vk} be the vector of the 
velocity components parallel to the line of sight at the spatial 
points Sj , determined by the recurrent sequence (see Paper 
II) 



;( Sj _i+As) = Wl- P+fvj-i , 3 = 2, 3, . . . , k .(1) 



Here 



/ = /(As) = exp(-|As|/0 



(2) 



is the correlation function and £ is a random number picked 
from a normal distribution with (£) = and Var(£) = 1. 

For a chosen 9 and velocity field distribution the objec- 
tive function C is calculated according to 



m 



- r(Ai) 



(3) 



where r(\i) is the simulated random intensity [eq.(19) in 
Paper II], the observed normalized intensity within 
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Figure 1. (a) - A Monte-Carlo simulation of the H+D Lyo profile. The thicker curve marks the blue wing of the HI line and the 
complex structure of the D absorption, (b) — The portion of the H+D profile marked in panel a is shown by open circles having error 
bars corresponding to S/N = 30 at the continuum level, while the solid curve is the result of the RMC minimization, (c) - The adopted 
one-dimensional velocity distribution (dotted curve) and the one determined by the RMC procedure (solid curve), (d) - Histograms p(u) 
for the adopted (dotted lines) and estimated (solid lines) projected velocity distributions (u = v/at)- 



the ith pixel of the line profile, at an estimate of the exper- 
imental error, and v = m — n the degree of freedom (m 
is the number of data points and n is the number of fitted 
physical parameters, n = 5 in our case). 

The minimization of \ 2 is the main goal of the proce- 
dure. Any parameter configuration with \ 2 ~ 1 has to be 
considered as a physically reasonable model for the interpre- 
tation of the observational data. 

The computational RMC procedure is split into two 
steps. At first, random values for the physical parameters 
are chosen, i.e. the vector 6. Secondly, an optimal velocity 
field configuration is estimated for these parameters. The 
current \ 2 value is computed and if it is larger than Xv,a 
(which is the value of \ 2 f° r a given credible probability 
P a — 1 — a) the whole procedure is repeated. In principle, 
to choose the vector 6 in the physical parameter subspace we 
could use any of the common optimization methods (like e.g. 
nonlinear simplex) but we apply the same stochastic RMC 
procedure for both stages of the computational scheme. In 
detail it is described as follows. 

1) A simulation box in the parameter subspace is specified 
by fixing the parameter boundaries. 

2) 8 is chosen arbitrarily in the simulation box. 

Next, the points 3-5 describe the evaluation of optimal ve- 
locity field configurations for a given set of physical param- 
eters. 



3) A random realization of the velocity field is generated, 
and the corresponding value of \ 2 is calculated. 

4) Subsequently an element Vj is chosen at random, and a 
random change is given to this element in accord with (Q). 
If j happens to be equal to 1, we put Vi = £o"t. 

5) The new value of C is calculated. If AC < 0, i.e. if 
Xnew < Xoidi the new 'trial' configuration of the velocity 
field is accepted. If AC > 0, the trial configuration is ac- 
cepted with a probability Pac = exp(— ^AC), i.e. we take 
a random number £ uniformly distributed in (0,1), and if 
( < Pac, the new configuration is accepted. If £ > Pac, 
the change is rejected and the configuration of the velocity 
field is returned back to its previous state. Then we repeat 
the procedure from step (4) up to /3i times to find the best 
solution (i.e. the smallest \ 2 value) for the 8 under consid- 
eration. Pi is chosen such that each Vj has been changed 
on average so often that Var(vj) ~ of . - The procedure is 
stopped when \ 2 X 2 ,a ■ 

6) An element of 8 is chosen at random, a random change 
A8i is given to it, and the procedure starting from step (3) 
is repeated (up to times) until Xmin wn l be less than 
an expected value xi,a ■ If after repetitions the condition 
Xmin < xt,a ls n °t satisfied the procedure stops. p2 is chosen 
such that each 8i has been changed on average no less than 
30 times to get a representative statistics. 
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Note that the convergence of the scheme depends on the 
size of the simulation box. The necessary size is unknown a 
priori. If it is too large, the computing time increases con- 
siderably, and if it is too small, the required x 2 value may 
not be reached at all. Therefore the box size must be chosen 
with some care and should be adjusted to the experimental 
data. 



2.2 Validity of the RMC simulations 

In order to test the validity of our RMC-computer code, we 
performed several tests in which the 'observational' H+D 
Lya profile was generated for a given set of the parameters 
§ = { N m = 2.5 xlO 17 cm~ 2 , D/H = 1.5 x 1(T 4 , T kin = 1250 
K, a t /vf h = 2.0, L/l = 5.0 }, and for one random realization 
of the one- dimensional velocity distribution. The chosen pa- 
rameters are typical for the warm diffuse interstellar gas (cf. 
Linsky et al. 1995) except for the D/H ratio which was taken 
approximately 10 times larger than the mean ISM value in 
order to make the deuterium absorption being more pro- 
nounced in the simulated spectrum. 

The test was conducted with the aim of evaluating 
whether the code is capable of finding the correct 6 and 
the correct structure of the stochastic velocity field. At first 
we considered a very soft constraint on the possible shape 
of the v(s) distribution. Namely, only the blue wing of the 
H+D Lya blend was used in the minimization procedure 
(this may correspond to the case when the red wing of the 
hydrogen profile is distorted by the Lya forest lines). We 
will see later (Section 3) that the uncertainty range of the 
D/H ratio determination found under this soft constraint 
may be considerably reduced by accounting for the shape of 
the whole H+D Lya absorption. 

Fig. la shows our simulated spectrum of H+D Lya. The 
example exhibits a complex intensity structure in DI. This 
complexity is caused by the finite correlation length effect 
only and, hence, it should not be interpreted as being caused 
by several (at least three) independent cloudlets with differ- 
ent radial velocities and different physical parameters. The 
velocity field distribution yielding this spectrum is shown in 
Fig. lc by the dotted curve. 

In order to illustrate the effects the correlated struc- 
ture of the velocity field has on the line forming process, 
Fig. 2 shows the local absorption coefficient k\(s) for HI 
and DI Lya for three fixed wavelength A ~ 1215.28, 1215.32 
and 1215.36 A which correspond to the two minima and the 
maximum between them in the DI profile shown in Fig. la. 

The correlations have the effect that for a given value 
of A the absorption coefficient may vary strongly along the 
line of sight. At the wavelengths chosen this effect is most 
pronounced for DI. This is in strong contrast to the micro- 
turbulent model (having the same a t and T kin ) in which the 
whole region along the line of sight would evenly contribute 
to the absorption at a given A. This is clearly seen from the 
hydrogen absorption coefficients which at the chosen wave- 
lengths, lie in the range of the damping wing where both 
micro- and mesoturbulent models have the same solutions 
(see Paper II). 

Let us turn back to Fig. 1. The thicker solid curve in 
panel a marks the portion of the H+D Lya profile used in 
the minimization procedure. To run the RMC-code we sam- 
pled the intensity at evenly spaced intervals AA and added 
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Figure 2. The contribution from different volume elements 
to the HI and DI local absorption coefficients for three fixed 
wavelengths A = 1215.279 A (al,a2), 1215.316 A (bl,b2), and 
1215.363 A (cl,c2), which correspond to the two minima and the 
maximum inbetween in the DI profile shown in Fig. la. 



to the /(A;) the 'experimental' errors which correspond to a 
signal-to- noise ratio S/N = 30 at the continuum level. The 
corresponding intensities and their uncertainties are shown 
in Fig. lb by the open circles and error bars. We consider 
these data as being 'observational' intensities and apply our 
RMC-code to estimate the set of n — 5 parameters { 9*j } 
and the corresponding velocity field configuration { v* }. 

The estimated model parameters are accepted if Xmin 
is found to rank in the xl,a distribution, i.e. Xmin — 1- The 
smooth curve in Fig. lb shows one of the possible RMC 

.2 



(cf. xi 



1.459 for v 



30 and a = 0.05). It corre- 



sponds to the model parameters N^j = 2.54 x 11T' cm 
(D/H)* = 1.44 x 10" 4 , n m = 1216 K, {o t /vf h Y = 2.02, 
(L/l)* — 5.06, and the ^-configuration shown in Fig. lc 
by the solid curve. It is worthwhile to emphasize once more 
that the derived ^-configuration is not unique : depend- 
ing on the initial condition (configuration of the velocity 
field) we may obtain several patterns for the velocity field 
distribution that satisfy our acceptance condition (each of 
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Figure 3. Confidence region for the plane 'log N(HI) - log D/H' 
for the fitted profile from Fig. lb when the other parameters T k i n , 
at/v^ h , and L/l are fixed but the configuration of the velocity field 
is free. The adopted parameters and those found by the RMC pro- 
cedure are labeled by the filled circle and diamond, respectively. 

them has its mirror counterpart which is also acceptable, of 
course). But the projected velocity distribution p(u) has ap- 
proximately the same shape for all possible RMC solutions. 
An example of p(u) for the adopted and estimated velocity 
distributions is shown in Fig. Id by the dotted and solid his- 
tograms, respectively [here u means the normalized velocity 
u(s) = v(s)/a t ]. 

The RMC procedure may result in the velocity field 
configuration becoming stuck in a local minimum of the C- 
function for a very long period ('meta-stable' state). Such a 
minimum may not solve the required problem at all. This dif- 
ficulty may be circumvented by making a few starts from dif- 
ferent random initial configurations. Numerical experiments 
have shown that for a given point 6 the Xmin values may 
scatter by about 30 per cent around the mean due to the 
probabilistic character of the solution for the velocity field. 
This fact should be taken into account when the error surface 
is calculated in order to determine the confidence regions for 
the estimated parameters { 8* }. To get meaningful confi- 
dence limits for fitted parameters the error surface must be 
calculated several times with different starting points. Hav- 
ing calculated a set of error surfaces, the best one may be 
obtained by simply picking the smallest \ 2 values for each 
point in the computational grid. 

Two examples of the confidence regions for the 'log Nm 
- log D/H' and 'log D/H - Tkin' projections are presented 
in Fig. 3 and 4, respectively. The results were obtained from 
a set of six computed error surfaces for each case. Shown 
are two confidence 'ellipses' - the regions constrained by 
the x 2 values corresponding to 68.3 % and 95.4 % C.L. (la 
and 2a regions for a normal distribution of residues, respec- 
tively). The RMC solutions are labeled by the filled dia- 
monds whereas the filled circles mark the adopted quantities 
used to produce the 'observed' spectrum. A good agreement 
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Figure 4. Confidence region for the plane 'log D/H - T kin ' for 
the fitted profile from Fig. lb when the other parameters JVhi, 
at /v?L , and L/l are fixed but the configuration of the velocity field 
is free. The adopted parameters and those found by the RMC pro- 
cedure are labeled by the filled circle and diamond, respectively. 

within statistical errors between these two points and be- 
tween the two projected velocity field distributions can be 
considered as a justification of our computational procedure. 

It is worth emphasizing that the elongated and declined 
shape of the confidence regions in Fig. 3 shows that the 
D/H value is strongly anti-correlated with TVhi. Therefore 
the D/H estimates may be badly dependent on the assumed 
hydrogen column density. 

The correlation between T kin and D/H is less pro- 
nounced. However, Fig. 4 shows that the kinetic temper- 
ature is not well determined in general. It seems reasonable 
to assume that the narrowest subcomponent of a complex 
absorption spectrum can show only an upper limit for Tkin 
since the correlated structure of the velocity field hampers 
the exact measurement of T kin from the apparent b-values 
(see Paper II, for details). 



3 APPLICATION TO QSO DATA 

The example from Subsection 2.2 demonstrates that com- 
plexity (and/or asymmetry) of an observed absorption spec- 
trum is not necessary an indication of a clumpy structure 
of the absorbing material but may be solely caused by the 
correlations in the velocity field. Therefore one may try to 
explain the observed spectra starting with a minimum num- 
ber of ad hoc assumptions concerning the structure of the 
absorber. In particular, one may consider a cloud model with 
a homogeneous density and temperature, assuming that the 
observed H+D absorption from QSO spectra arises in the 
outer regions (halos) of the putative intervening galaxies (cf. 
Tytler el al. 1996). If this is the case, then Nm ranges from 
~ 10 17 to ~ 10 18 cm~ 2 . This in turn means that the geo- 
metrical size of the absorbing material (parameter L in our 
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model) may be very large, L ~ 10 kpc. If there are many ran- 
dom velocity elements along the line of sight (L/l 3> 1), the 
theoretical expectation value should be a good approx- 
imation to the observations (Paper I) . If, however, the ratio 
L/l is not too large, deviations of from the observed 
profile (different kinds of asymmetry) are to be expected. 

Modern observations of the H+D absorption at high 
redshift revealed two limiting D/H ratios ~2x 1CF 4 (e.g. 
Carswcll et al. 1996), and ~ 2 x 10" 5 (e.g. Tytler et al. 
1996). The primordial D/H ratio is shown to be very sen- 
sitive to the physical conditions that existed at the time of 
the Big Bang (e.g. Walker et al. 1991). D/H is expected to 
decrease with cosmic time due to processing of gas by stars. 
The theoretical uncertainties in the calculations of the pri- 
mordial hydrogen isotopic ratio are rather small, the D/H 
value being predicted for a given baryon-photon ratio, ij, 
with an accuracy of ~ 15 % (see Sarkar 1996). Therefore, 
accurate measurements of the primordial D/H ratio would 
provide important constraints on the cosmological theories. 
Since these measurements are of fundamental importance 
to cosmology, the validity of the revealed two limiting D 
abundances and the reliability of the quoted errors should 
be thoroughly investigated. 

In the following subsection we will show that accounting 
of the velocity field structure of the large scale gas flows may 
have a major impact to the D/H measurements. 

3.1 The z a = 2.504 system in Q 1009+2956 

The RMC procedure was applied to the D absorption ob- 
served by B&T at z a = 2.504 in front of the quasar Q 
1009+2956. 

Using a iwo-cloud microturbulent model and applying 
a Voigt profile fitting procedure, B&T derived the following 
parameters for this system : total hydrogen column density 
N m = 2.9±g;| x 10 17 cm" 2 and D/H = 3.0±g;| x 10~ 5 ; indi- 
vidual T kin = 2.I+0.1 x 10 4 K and 2.4j£ 7 . x 10 4 K for the 
blue and red subcomponents, respectively, the correspond- 
ing turbulent velocities b tU rb = 3.2±0.4 km s _1 and 2.3±1.4 
km s _1 , and a difference in the radial velocity (derived from 
metal lines) of 11 km s~ . Including corrections for possible 
weak HI absorption at the position of D (estimated as an 
expected value in Monte Carlo simulations), B&T reduced 
the measured D/H ratio to log (D/H) = -4.60 ± 0.08 ± 0.06 
(i.e. D/H ~ 2.5 x 10~ 5 ), where the la random error followed 
by the systematic error from fitting the continuum level are 
given. However, for a given line of sight the proper correc- 
tion is never known, and therefore we will not reduce D/H 
values in our analysis. This should be borne in mind when 
discussing the true D abundance. 

We start with the calculation of a template spectrum 
using the parameters of B&T listed in their Table 1. To sim- 
ulate real data, the experimental uncertainties are added to 
the computed intensities which were sampled in equidistant 
bins as shown in Fig. 5(al,bl,cl) by dots and correspond- 
ing error bars. For the redshift the mean value of 2.504 was 
adopted and a one-component mesoturbulent model with a 
homogeneous density and temperature was used to construct 
the objective function. 

Adequate profile fits (Xmi« P er degree of freedom <, 1) 
for three different sets of parameters are shown in pan- 
els (al,bl,cl) of Fig. 5 by solid curves, whereas the in- 



dividual HI and DI Lya profiles are shown by dashed and 
dotted curves, respectively. The estimated parameters and 
corresponding xLn values are also listed in these panels 
for each solution. The adjusted t((a;)-configurations (pan- 
els a2,b2,c2) and their projected distributions p(u) (panels 
a3,b3,c3) show pronounced deviations from a single Gaus- 
sian. Such distributions may, in principle, have very long 
positive 'tails' yielding unobservable D absorption which is 
hidden by the core of the hydrogen line. However, including 
in the objective function both the blue and the red wing of 
the HI line, allows to constrain significantly the set of the 
u(a')-configurations. As a result the uncertainty for the D/H 
measurements is considerably reduced. 

In accord with our calculations, the restored profile 
of the DI line in the z a = 2.504 system is almost sym- 
metrical. A slight asymmetry of the red DI Lya wing is 
caused by the weak distortion of the p(u) distribution to- 
ward positive u-values. In this case the width of the DI 
line may be formally characterized by the 6-parameter : 
6 D i = vf h \\ + 2(a t /vf h ) 2 ] L/2 . The inspection of Fig. 5 
(al,bl,cl) reveals rather small bui value (~ 15 km s _1 ) 
as compared with a t ~ 25 km s _1 obtained by the RMC 
procedure. This is not a surprise, however, because random 
realization of the correlated velocity field may have p(u) 
which deviates substantially from the ensemble average dis- 
tribution (assumed to be Gaussian in this case) leading to a 
'contraction' of this distribution (see Paper II). 

One may further be interested in the question whether 
such weak asymmetry of p(u) is essential for the analysis of 
the H+D Lya profile, or whether a simpler one component 
microturbulent model could be used to study the z a — 2.504 
absorber. The answer is that it is essential indeed, and that 
the simplified model fails to reproduce the profile in ques- 
tion. The smallest Xmin value obtained for this model is only 
2.4 which rejects the one-component microturbulent solu- 
tion. This example clearly demonstrates once more that the 
analysis of absorption spectra may be very sensitive to the 
model assumptions concerning the velocity field structure. 

The essential difference between the results of B&T and 
ours lies in the estimation of the hydrodynamical velocities 
in the z a = 2.504 absorbing region. The two-component mi- 
croturbulent model yielded bf lue — &P ed = 14 km s _1 and 
bt U rb — 3 km s _1 which is about 10 times smaller than \f2a t 
estimated by the RMC procedure. For the additional pa- 
rameter L/l of the mesoturbulent model, we found a value 
of ~ 5, indicating that the space averaging along the line of 
sight corresponds to an incomplete statistical sample and, 
hence, significant deviations of the observed intensities from 
their expectation values may occur. 

As for the other parameters (Nm, D/H, and !«„), the 
RMC method yields in this particular case values which are 
similar to those of B&T. The total hydrogen column den- 
sities and the D/H ratios obtained are shown in Fig. 6 by 
the filled diamond for the best RMC solution (panel bl in 
Fig. 5) and by the filled circle for the B&T model (the er- 
ror bars shown correspond to la in accord with the B&T 
data from Table 1). The filled triangle and square mark the 
(Nm, D/H) pairs found separately for the blue and red com- 
ponents. The D abundance for the red component was as- 
sumed by B&T as (D/H) red = (D/H) w „ e = (D/H) toto , and, 
thus, there was no real measurements of the D/H value for 
the red component in their model. This is marked by arrows 
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Figure 5. (al, bl, cl) — A template H+D Lya profile (dots with error bars) representing the normalized intensities and their 
uncertainties in accord with the B&T data. The solid curves show the results of the RMC minimization, whereas the dotted and dashed 
curves are the separate profiles of DI and HI, respectively. Also shown are the best fitting parameters and the \ 2 values obtained for 
each case. (a2, b2, c2) - The corresponding individual realizations of the one-dimensional velocity distributions u(x) in units of at- (a3, 
b3, c3) - The histograms are the projected velocity distributions p(u). 



instead of error bars at the filled square in Fig. 6. The pro- 
jection of the 5-dimensional error surface onto the 'log Nm - 
log D/H' plane was calculated by the RMC procedure using 
the A(x 2 ) method. Drawn are 68 % and 95 % confidence 
regions computed under the condition that the parameters 
Thin, o-t/vth, and L/l are fixed, but the velocity field config- 
uration is free. As seen, the B&T value for the total A^hi is 
located closely to the RMC confidence regions, but the D/H 
ratio is systematically lower in the microturbulent model. 

Fig. 7 shows 68 % and 95 % confidence regions for 
the two parameters Thin and D/H computed under the as- 
sumption that JV H i, crt/vth, and L/l are fixed, but the u(x)- 
configuration is free. The kinetic temperature is found in the 



range from ~ 15000 K to ~ 22000 K, whereas D/H lies be- 
tween ~ 3.3 x 10~ 5 and ~ 4.5 x 10~ 5 . In this figure, the filled 
symbols label the RMC and B&T best fitting parameters in 
the same sense as in Fig. 6. Note the marginal position of 
the RMC solution (the filled diamond) which shows that 
there are points in the 'D/H - Tki„' plane with \ 2 values 
less than Xmin = 0.485 adopted for the best RMC solution 
at the first stage of calculations. Fig. 7 shows that also in 
the 'D/H - Thin plane the D/H uncertainty region is shifted 
toward higher D/H values as compared with the B&T data. 

Above it was noted that the D/H ratio does sensitively 
depend on the Am value chosen. The measured D/H value 
and its confidence interval are both also affected by the 



© 0000 RAS, MNRAS 000, 000-000 



8 S. A. Levshakov et al. 



-4.3 



-4.4 - 



S3 

Q 



OA 

O -4 5 
— 1 



-4.6 



▲ 



1 i 1 r 

68.3 % C.L. 

95.4 % C.L. 

RMC - best solution 
B&T - red component 
B&T - blue component 
B&T - total 




16.8 



17.0 



17.2 



17.4 



Log N(HI), cm 



Figure 6. Confidence region for the plane 'log N(HI) - log D/H' 
for the fitted profile from Fig. 5(bl) when the other parameters 
Tfeim a t/v^ h , and L/l are fixed but the configuration of the veloc- 
ity field is free. The RMC solution is labeled by the filled diamond. 
The filled square and triangle represent the B&T parameters for 
the red and blue components of their two-component microturbu- 
lent model, and the filled circle marks the combined solution. The 
error bars correspond to la in accord with B&T, but the arrows 
reflect the fact that the D abundance was estimated through the 
blue component mainly and (D/H) re£ j = (D/H)biue was assumed 
by B&T. 



30000 



25000 



20000 



15000 



I 1 I I I 1 I 


!>>>!>>>!>>>!>> 

68.3 % C.L. 

95.4% C.L. 

♦ RMC - best solution 
■ B&T - red component - 
A B&T - blue component - 


- — 1 


1 — 


- i i 

1 , , , 1 , 


, I , , , iT-rTT'r, I , , 



2.4 



2.8 



3.2 



3.6 



4.0 



4.4 



(D/H) x 10 

Figure 7. Confidence region for the plane 'D/H - Tfcj„' for 
the fitted profile from Fig. 5(bl) when the other parameters Nm, 
at I v^ h , and L/l are fixed but the configuration of the velocity field 
is free. The symbols used have the same meaning as in Fig. 6. 



choice of the velocity field configuration. The accuracy of 
the RMC solution may be increased if some additional in- 
formation about the velocity field structure is available for 
the analysis. For instance, the profile shapes of higher order 
Lyman lines may yield additional useful constraints to the 
allowable p(u) distribution. For this case the size of the D/H 
uncertainty interval may be narrowed down. This is demon- 
strated in Fig. 8 where the results of two calculations (for- 
mally called 'high' and 'low' D/H) are shown. The 68 % and 
95 % confidence regions were computed for the fixed Tki n , 
a t /v th , and L/l listed in panels cl and al (Fig. 5), respec- 
tively, and for the fixed configurations of u(x) shown in pan- 
els c2 and a2 (Fig. 5). Provided the constraints mentioned 
above are available and, thus a distinction between differ- 
ent p(u) distributions compatible with the measured H+D 
profile (Fig. 5) is possible, the accuracy of the D/H mea- 
surements improves essentially becoming comparable with 
the accuracy of the SBBN calculations (Sarkar 1996). 

We now discuss the observed D abundance and its in- 
ferred primordial value which can be compared with the 
SBBN model predictions. Within the framework of the 
SBBN, the primordial abundances of the light elements (D, 
3 He, 4 He, 7 Li) depend only on the baryon to photon ratio 
f) = nt/n 7 , where rib is the current cosmological number 
density of baryons and n 7 is the number density of the cos- 
mic background photons. The predicted relations between 
the abundances of D, 4 He, and 7 Li and r\ are illustrated in 
Fig. 9. The thick solid curves plotted in this figure are com- 
puted by using the analytical formulae from Sarkar (1996), 
while dashed curves depict their uncertainty intervals. The 
primordial mass fraction of 4 He, Y p = 0.243T 0.003, was de- 
rived by Izotov et al. (1997), and Bonifacio & Molaro (1997) 
give the following value for the lithium abundance in halo 
population II stars : log ( 7 Li/H)^ = -9.762 ± 0.012 ± 0.05 
(here the subscript p denotes primordial). The uncertainty 
regions for the measured values of Y p ( 4 He) and ( 7 Li/H)p : 
and the implied SBBN limits on r\ are shown by the solid- 
line rectangles in Fig. 9(b,c). The upper panel a in this figure 
gives our results : the solid-line rectangle is defined by the 
95 % confidence region (shown in Fig. 6) for an arbitrary 
configuration of the velocity field, whereas two dotted-line 
rectangles correspond to 68 % confidence regions computed 
for the two fixed velocity fields leading to the 'high D/H' 
and 'low D/H' solutions shown in Fig. 8. The shaded re- 
gion in Fig. 9 represents the range of r\ compatible with the 
measured abundances of all three nuclei. Fig. 9 shows that 
concordance of D/H with Y p ( 4 He) and ( 7 Li/H)p measure- 
ments is achieved if 3.9 x 10~ 10 £ rj £ 5.3 x 10~ 10 . This 
constraint on rj is consistent with results of other groups; 
for example, Hata et al. (1995) quote r) = 4.4±°;§ x 10~ 10 
(la) [see also Sarkar (1996) for a discussion of the inter- 
comparisons of the allowed ^-intervals] . The center of the 
shaded ^-interval in Fig. 9 corresponds to O^/i 2 = 0.017, 
where Q b is the fraction of the critical density contributed 
by baryons and h the Hubble constant in units 100 km s _1 
Mpc -1 . The uncertainties in this estimation can alter fib by 
up to ~ 20 %. 
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4 CONCLUSIONS 

The origin of QSO absorption-line systems with HI column 
densities ~ 10 17 - 10 18 cm (which are particularly use- 
ful for deuterium observations) is usually associated with 
extended gas halos of foreground galaxies. According to re- 
cent high resolution observations, the line profiles in these 
systems are often wider than expected from purely thermal 
broadening, and, therefore, their widths reflect both ther- 
mal and turbulent broadening caused by any kind of bulk 
motions in the halos. In turbulent media the determination 
of HI and DI column densities depends on the model (micro- 
or mesoturbulent) is chosen for the profile analysis. The mi- 
croturbulent model yields reasonable results if a t <C Vth , or 
L/l 3> 1, otherwise the mesoturbulent one is more adequate. 

Direct observations of the Lya emission of distant 
galaxies [z > 2) show that their giant halos (R ~ 100 kpc) 
have complex morphologies and kinematics (Rottgering et 
al. 1996). The observations reveal as well extended HI ab- 
sorption systems (with projected sizes up to ~ 50 kpc) seen 
against the Lya emission (van Ojik et al. 1997). 

Using the measured Doppler parameters for the absorb- 



ing gas with Nm 



If) 1 



10 1 



cm (see Table 3 in 



van Ojik et al.) and the kinetic temperature T^ in = 10 4 K 
adopted by van Ojik et al., we can estimate the range for the 
velocity dispersion of bulk motions : 24 km s _1 £ a t S5 55 
km s _1 . It follows that at/vth > 1 within the HI absorp- 
tion line gas embedded in the galactic halos. Both findings 
- high absolute value of at and high at/vth ratio, - make 
the mesoturbulent approach to appear to be more realistic 
one. Besides, our present RMC calculations yielded for the 
z a — 2.504 absorption system a t — 25 km s _1 which lies just 
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Figure 9. Comparison of predicted primordial abundances (de- 
noted by the subscript p) with observational bounds. The the- 
oretical SBBN D, 4 He, and 7 Li yields (solid curves) and their 
uncertainties (dashed curves) as function of r] (the baryon to pho- 
ton number ratio) are computed by using the parameterization of 
Sarkar (1996). The abundances of D and 7 Li are number ratios, 
whereas Y p is the mass fraction of 4 He. The hight of the solid-line 
rectangles in panels b and c give the bounds from recent measure- 
ments of extra-galactic Y p (Izotov et al. 1997) and 7 Li in popu- 
lation II (halo) stars (Bonifacio & Molaro 1997). The horizontal 
widths give the range of r\ compatible with the measurements. In 
panel a the solid-line rectangle corresponds to the la uncertainty 
region for D/H towards Q1009+2956, as calculated by the RMC 
procedure for an arbitrary velocity field configuration (see Fig. 6), 
whereas the two dotted-line rectangles are defined by the 68 % 
confidence regions shown in Fig. 8 (the velocity field structure 
is fixed). The shaded region is the window for r/ common to all 
results. 
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within the observed range, whereas at — 2 km s found by 
B&T is evidently too low. 

The main conclusion of our work is that the mesoturbu- 
lent approach together with the RMC computational scheme 
allows to solve the inverse problem for the H+D Lya absorp- 
tion and to restore the physical parameters of the absorbing 
gas as well as the projection of the velocity field distribution. 

The proposed computational procedure enables us to 
draw confidence regions for different pairs of the adopted 
parameters. From the analysis of the synthetic H+D Lya 
spectrum it was found that, in general, the measured D/H 
and Nm values are anti-correlated. 

The study of the template H+D Lya profile (which 
reproduces the original spectrum of B&T) yields D/H = 
(3.75 ± 0.85) x 1CT 5 (2a). We therefore adopt the value of 
4.6 x 10~ 5 as a conservative upper limit on the primordial 
abundance of D relative to hydrogen. 

We conclude that the discordance of D/H with the 4 He 
and 7 Li primordial abundances noted by B&T is a conse- 
quence of the use of the microturbulent model. Within the 
framework of the generalized model one finds good agree- 
ment between the measurements mentioned above and the 
SBBN predictions. 
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